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Abstract 

Stability  is  proven  for  two  second  order,  two  step  methods  for  uncoupling  a  system  of 
two  evolution  equations  with  exactly  skew  symmetric  coupling:  the  Crank-Nicolson 
Leap  Frog  (CNLF)  combination  and  the  BDF2-AB2  combination.  The  form  of  the 
coupling  studied  arises  in  spatial  discretizations  of  the  Stokes-Darcy  problem.  For 
CNLF  we  prove  stability  for  the  coupled  system  under  the  time  step  condition 
suggested  by  linear  stability  theory  for  the  Leap-Frog  scheme.  This  seems  to  be  a 
first  proof  of  a  widely  believed  result.  For  BDF2-AB2  we  prove  stability  under  a 
condition  that  is  better  than  the  one  suggested  by  linear  stability  theory  for  the 
individual  methods.  This  report  is  an  expended  version  of  the  one  submitted  for 
publication. 
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1  Introduction 


This  is  an  expanded  version,  containing  supplementary  material,  of 

a  report  with  the  same  title. 

In  this  report  we  prove  stability  of  two,  second  order  IMEX  methods  for 
uncoupling  two  evolution  equations  with  exactly  skew  symmetric  coupling: 
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du 

— — h  A\u  +  C(j)  =  f{t),  for  t  >  0  and  n(0)  =  uq 

LLL 

+  A2(j)  —  CTu  =  g(t ),  for  t  >  0  and  0(0)  =  0o- 

LLL 


This  problem  occurs,  for  example,  after  spatial  discretization  of  the  evolution¬ 
ary  Stokes-Darcy  problem,  e.g.,  [16,12,18,17].  Here 


u  :  [0,  oo)  — >  M.N ,  0  :  [0,  oo)  — >  RM, 


and  f,g,uo,(fto  and  the  matrices  Ai/2,C  have  compatible  dimensions  (and  in 
particular  C  is  N  x  M).  Note  especially  the  exactly  skew  symmetric  coupling 
linking  the  two  equations.  We  assume  that  the  Ai  are  SPD.  Our  analysis 
extends  to  the  case  of  Ai  positive  real  or  even  nonlinear  with  (A(v),v)  > 
Const. |v|J.  With  superscript  denoting  the  time  step  number,  the  first  method 
is  CNLF,  the  combination  of  Crank-Nicolson  and  Leap  Frog  given  by:  for 
n  >  2 
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(CNLF) 


Since  the  stability  region  of  LF  is  the  interval  —1  <  Im(z)  <  +1,  from  the 
scalar  case  we  expect  a  stability  restriction  of  the  form  A tyj Amax(C,TC)  <  1. 
Interestingly,  it  seems  that  sufficiency  in  the  non-commutative  case  is  not  yet 
proven  Verwer  [22],  remark  3.1,  page  6.  We  prove  in  Section  2  that  CNLF  is 
indeed  stable  under  (1),  exactly  the  condition  suggested  by  the  linear  stability 
theory. 

For  vectors  of  the  same  length,  denote  the  usual  euclidean  inner  product  and 
norm  by  (u,v)  :=  uTv  ,  |0|2  :=  (0,0).  We  denote  the  weighted  norms  by 


mIax  :=  «TAi-u,  and  |0|^2  :=  0TH20. 


Theorem  1  (Stability  of  CNLF)  Consider  CNLF.  Suppose  the  time  step 
restriction  holds: 


At\J\maACTC)  <  a  <  1,  for  some  a  <  1. 
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Then  for  any  n  >  2 
1  —  ol  r 
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Next  we  establish  the  stability  of  BDF2  with  explicit  AB2  coupling 
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(BDF2-AB2) 


The  stability  region  of  AB2  suggests  that  this  combination  is  strictly  worse 
than  CNLF.  However,  we  prove  that  the  combination  inherits  enough  stability 
from  BDF2  to  be  stable  under  a  time  step  condition  that  in  many  cases  is 
better  than  the  one  for  CNLF. 

Theorem  2  (Stability  of  BDF2-AB2)  Consider  BDF2-AB2.  Suppose  that 
the  time  step  restriction  holds 

At  max{Amax(Ab1CCT),  AmaAAf  1CTC)}  <  a  <  1,  for  some  a  >  0,  (2) 

then  BDF2-AB2  is  stable: 

\un\2  +  |0n|2  <  Cfinitial  data,  forcing  terms),  for  any  n  >  2. 

More  precisely,  for  all  n  >  1,  we  have  that 

t(|f”+1r+|f“+1|2)  +  i(|2u"+W|2+  \2<j>n+l—cf>n\‘I)  +  (K(+1+'H'+1 

<  i(l«‘|2  +  iy  I2)  +  i (12m1  -  «T  +  |2 01  -  /I2) 

+Aty _ 1 _ ( +  l3wl2  i, 
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where  we  have  denoted 
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Note  that  (2)  implies  that  A\  —  AtCTC,  A 2  —  A tCCT  are  SPD. 

Both  methods  use  3  levels;  approximations  are  needed  at  the  first  two  time 
steps  to  begin.  We  suppose  these  are  computed  to  appropriate  accuracy,  Ver- 
wer  [22], 

Because  the  problem  and  methods  are  linear,  stability  immediately  implies 
that  the  error  is  bounded  by  its  consistency  error. 


1.1  Connection  to  the  coupled  Stokes-Darcy  problem 


To  specify  the  motivating  problem  leading  to  the  system  of  evolution  equations 
considered,  let  two  domains  be  denoted  by  flf,  klp  and  lie  across  an  interface  / 
from  each  other.  The  fluid  velocity  and  porous  media  piezometric  head  (Darcy 
pressure)  satisfy 


ut  —  vAu  +  Vp  =  if(x,  t).  V  •  u  =  0,  in  Qf,  (3) 

S^t  -  V  •  (/CV0)  =  fp,  in  np, 

0)  =  0o,  hr  C.p  and  u(x,  0)  =  no,  in  Qf, 

4>(x,t)  =  0,  in  dklp\I  and  u(x,t )  =  0,  in  dflf\I, 

+  coupling  conditions  across  /. 

Let  hf/p  denote  the  indicated,  outward  pointing,  unit  normal  vector  on  I.  The 
coupling  conditions  are  conservation  of  mass  and  balance  of  forces  on  / 

u  ■  hf  +  Up  •  rip  =  0,  on  I  u  ■  hf  —  -ICVcJ)  ■  hp  =  0,  on  1, 
p  —  v  hf  ■  V«  •  hf  =  pg<f>  on  1. 

The  last  condition  needed  is  a  tangential  condition  on  the  fluid  region’s  velocity 
on  the  interface.  The  most  correct  condition  is  not  completely  understood 
(possibly  due  to  matching  a  pointwise  velocity  in  the  fluid  region  with  an 
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averaged  or  homogenized  velocity  in  the  porous  region).  We  take  the  Beavers- 
Joseph-Saffman  (-Jones)  interfacial  coupling 

^  a  ^  ^ 

—v  Vu  ■  rif  =  ==tt  ■  Tii  on  I  for  any  rt  tangent  vector  on  I. 

V  7~i  *  /C  *  Ti 

This  is  a  simplification  of  the  original  and  more  physically  realistic  Beavers- 
Joseph  conditions  (in  u  ■  %  which  is  replaced  by  (u  —  up)  •  %).  Here: 


(f>  =  Darcy  pressure  +  elevation  induced  pressure  =  piezometric  head 
q=  volume  discharge, 

Up  =  fluid  velocity  in  porous  media  region,  Qp, 
u  =  fluid  velocity  in  Stokes  region, 
if,  /p  =  body  forces  in  fluid  region  and  source  in  porous  region, 

/C  =  hydraulic  conductivity  tensor, 
v  =  kinematic  viscosity  of  fluid, 

S'o  =  specific  mass  storativity  coefficient, 

7]  =  volumetric  porosity, 
p  =  density, 

g  =  gravitational  acceleration  constant. 

We  shall  assume  that  the  boundary  conditions  are  simple  Dirichlet  conditions 
on  the  exterior  boundaries  (not  including  the  interface  I). 

We  denote  the  L2(I)  norm  by  jj  •  || /  and  the  L2(Qf/p)  norms  by  jj  •  \\f/p, 
respectively;  the  corresponding  inner  products  are  denoted  by  (•,  •)//??•  Define 


Xf  :  ={ve  (iJ1(D/))fl!  :  v  =  0  on  d Qf\I}, 
Xp  :  =  {V>  e  H\np)  :tj}  =  0  on  dQp\I}, 

Q  =  Ll(Slf). 


Define  the  bilinear  forms 


_ ^  r  qi 

af(u,v)  =  (Wu,  Vu)/  +  V  /  ^  (u  •  ?i)(v  •  ?i)ds, 

i  ^  I  \/  7~i  *  /C  • 

ap(0,^)  =  (A'V</>,  V^)p, 

Cj(u,  (p)  =  npg  J  (pu-fifds. 

A  (monolithic)  variational  formulation  of  the  coupled  problem  is  to  find  (u,  p,  (p)  : 
[0,  ex:)  — >  Xf  x  Qj  x  Xp  satisfying  the  given  initial  conditions  and,  for  all 
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(4) 


v  E  Xf,  q  E  Qf,  0  E  Xp 

(' ut i  v)f  +  af(u,  v )  -  (p,  V  •  v)f  +  a(v,  4>)  =  {if,  v)f, 

{q,  v  -u)f  =  o, 

So((fH,  ip)p  +  dp{(f>,  0)  -  Cj(u,  0)  =  ( fP ,  ij})p. 

Note  that,  setting  v  —  u,ip  —  4>  and  adding,  the  coupling  terms  exactly  cancel 
in  the  monolithic  sum  yielding  the  energy  estimate  for  the  coupled  system. 

To  discretize  the  Stokes-Darcy  problem  in  space  by  the  finite  element  method, 
we  select  finite  element  spaces 

velocity:  X}j  C  Xf,  Darcy  pressure:  Xp  C  Xp,  Stokes  pressure:  Qhf  C  Qf 

based  on  a  conforming  FEM  triangulation  with  maximum  triangle  diameter 
denoted  ”/i”.  No  mesh  compatibility  at  the  interface  /  between  the  FEM 
meshes  in  the  two  subdomains  is  assumed.  The  Stokes  velocity-pressure  FEM 
spaces  are  assumed  to  satisfy  the  usual  discrete  inf-sup  condition  for  stability 
of  the  discrete  pressure.  We  denote  the  discretely  divergence  free  velocities  by 

Vh  :=  Xhf  n  K  :  (qh,  V  •  vh) f  =  0,  for  all  qh  E  Qhf} 

The  semi-discrete  approximations  are  maps  ( Uh,Ph ,  <f>h)  '■  [0;  oo)  — >  Xf  xQj  x 

X^  satisfying  the  given  initial  conditions  and,  for  all  ty  E  Xf,  qh  E  Qhf,f>h  £ 
yh 
V 

{Uh,t,vh)f  +  af{uh,vh)  -  (ph,v  ■  vh)f  +  c/(uh,0h)  =  (f f,vh)f, 

(qh,  V  •  uh)f  =  0,  (5) 

*5*0 ^h)p  T  Ojp{4^hi  Iph)  Ci{Uh,  'tph)  ( fp i  0 h)p ■ 

Note  in  particular  the  exactly  skew  symmetric  coupling  between  the  Stokes 
and  the  Darcy  sub-problems.  If  the  velocity  is  restricted  to  the  discretely 
divergence  free  subspace  we  obtain  The  semi-discrete  approximations  are  maps 
(uh,  4>h)  :  [0,  oo)  — >  Vf1  x  X^  satisfying  the  given  initial  conditions  and,  for  all 
vh  E  Vfi,'^h  E  Vp 

(uh,uvh)f  +  af(uh,vh)  +cI{vh,(/)h)  =  (f f,vh)f, 

So {4*h,ti  ifthjp  T  Up(0/i,  0/i)  Cj{llh,  0/i)  ( fp ,  1ph)p- 

The  exactly  skew  symmetric  coupling  between  the  Stokes  and  the  Darcy  sub¬ 
problems  is  retained.  Picking  a  basis  for  the  FEM  spaces  in  the  above,  this 
leads  to  a  system: 

du 

M/—  +  A\u  +  C(f>  =  f(t),  for  t  >  0  and  w(0)  =  uq 
S{)Mp-^-  +  A2(p  —  CTu  —  gif),  for  t  >  0  and  0(0)  =  0o- 

LLL 
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Here  the  respective  FEM  mass  matrices  are  denoted  Mf/p.  These  are  often 
spectrally  equivalent  to  the  identity.  The  above  system  can  also  be  reduced 
to  the  one  studied  by  a  further  change  of  variable.  The  Stokes-Darcy  problem 
has  experienced  a  rapid  development  of  numerical  methods.  We  end  the  paper 
with  a  list  of  some  additional  papers  that,  while  not  relevant  to  the  precise 
problem  considered  herein,  are  on  numerical  methods  for  the  Stokes  -Darcy 
problem. 


1.2  Previous  work 


When  Ai  are  SPD,  IMEX  methods,  like  CNLF  and  BDF2-AB2  require  the  so¬ 
lution  of  two,  smaller  SPD  systems  per  time  step  (which  can  be  done  by  legacy 
codes  for  the  independent  sub-problems)  as  compared  to  one  larger,  nonsyrn- 
metric  system  for  monolithically  coupled  methods.  Given  this  potentially  large 
simplification,  it  is  not  surprising  that  IMEX  methods  (and  associated  parti¬ 
tioned  schemes)  have  been  used  extensively  in  the  computational  practice  of 
multi-domain,  multiphysics  applications.  The  theory  of  IMEX  methods  is  also 
developing;  see  [11,21,2]  and  [13]  for  early  papers  and  [1,10]  and  particularly 

[22]  and  the  book  [14]  for  recent  work.  CNLF  is  itself  a  classic  (e.g.  [15])  com¬ 
bination  of  methods  in  computational  fluid  dynamics  with  wide  practical  use, 
including  in  the  dynamic  core  of  the  NCAR  climate  model,  [19]. 

Partitioned  methods  are  often  motivated  by  available  codes  for  subproblems 
[20]  and  tend  to  be  application  specific.  Examples  of  partitioned  methods  in¬ 
clude  ones  designed  for  fluid-structure  interaction  [3,4,7],  Maxwell’s  equations 

[23]  and  atmosphere-ocean  coupling  [8,10,9].  The  block  system  we  study  arises 
in  evolutionary  groundwater-surface  water  coupling,  e.g.,  [6,5,12,16].  Mu  and 
Zhu  [17]  gave  the  first  (in  2010)  numerical  analysis  of  a  partitioned  method 
based  on  the  backward  Euler-forward  Euler  IMEX  scheme;  this  has  been  ex¬ 
tended  to,  so-called,  asynchronous  time  stepping  (different  time  steps  for  dif¬ 
ferent  system  components)  in  [18].  Our  work  herein  is  motivated  by  the  search 
for  partitioned  methods  for  the  Stokes-Darcy  problem  with  higher  accuracy 
and  better  stability. 


2  Proof  of  stability  of  CNLF 


This  section  gives  a  complete  proof  of  Theorem  1. 
Lemma  3  We  estimate 

(C'4>,u)  -  f|C0|2  +  h«|2  -  1|m  -  CV|2 
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and,  if  Ai  are  SPD 


\u\  <  101  <  Amf(A2)|0U2, 

\c<f>\  <  yJXmzxicrcM. 

Thus 

\(C<t>,u)\  <  l^\mUCTCM2+y\max(CTC)\u\2. 

Proof.  The  first  claim  is  the  polarization  identity.  The  second  inequality  is 
elementary  while  the  fourth  follows  by  inserting  the  third  into  the  first.  For 
the  third,  we  have 

|CV>|  =  { 2  =  (CtC4>A)112  <  Kil(CTC)\^\. 


The  first  of  three  main  steps  in  the  proof  of  Theorem  1  is  to  take  the  inner 
product  of  CNLF  with  un+1  +  un~ 1  and  0n+1  +  0n_1  and  add: 
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The  second  step  is  to  rearrange  the  coupling  terms  as  an  exact  difference 
between  two  time  levels:  Coupling  =  (C(f>n,un+1  —  un~l)  —  (CTun,(pn+1  — 
077-1)  =  Cn+ 1/2  -  A”-1/2,  where 


An+1/2  :  =  (Cfn,un+1)  -  (C(j)n+1,un), 
Cn~1/2  :  =  (A0n_1,un)  -  (A0n,  un_1). 


The  third  step  is  to  add  and  subtract  \un\2  +  \  <f>n\ 2  to  the  control  the  energy 
at  level  tn: 
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=  (/",  un+i  +  tt""1)  +  (0\  un+1  +  un~l)  =  RHS. 
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Using  Lemma  3  we  treat  RHS  in  a  standard  way: 


RHS  <  |r|Am‘„/2(Ai)|«“+1  +«"-1U,  +  l9”|Ai/2(A2)|r+1  +  F-'\m 

<  (KiiAMr+Ki^m2)  +  \(w,+i+u"-'\x+\r+1+r-l\%)- 


Thus,  define  the  system  energy 


j^n+1/2 


1 

2  L 


u 


n+1 1 2 


+  ir+r  +  iuT  +  \<t> 


in\2 


+  A  tCn+1/2. 


Collecting  terms  we  obtain 

En+ 1/2  -  En+1/2  +  At  (\un+1  +  un~l\2M  +  | cj)n+1  +  0n_1|^2) 

<At(A-Jn(Ai)in2  +  A-L(dl2)|^f). 


Obviously,  En+ x/2  —  A”-1/2  +  {positive_terms}  <  RHS  immediately  implies 
stability  provided  only  that  En+1/2  >  0  for  every  n.  We  have  (using  Lemma 
3  to  bound  the  coupling  terms) 


7771+1/2  >  _ 

-  2  L 


U 


n+1\2  +  \(j)n+1\2  +  \un\2  +  \(f)n\2 


- 7T\I  Amax(C,TC) 


U 


n+1'2  +  \un\2  +  \(j)n+1\2  +  \(j) 


n\2 


This  is  positive  (completing  the  proof)  provided 


At\J Amax(ATA)  <  1. 


3  Proof  of  stability  of  BDF2-AB2 


We  proceed  to  prove  Theorem  2.  Take  the  inner  product  of  BDF2-AB2  with 
un+l,  0n+1,  respectively,  and  add.  There  are  two  keys  to  the  proof  of  stability. 
The  first  key  is  the  treatment  of  the  BDF2  term.  Apply  the  identity 


(2a  —  6)2j 

4  . 


b2  (2b  -c)2^ 
J 4  4 


+ 


2b  +  c)s 


4  b  +  c)a 
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with  a  =  un+1,b  =  un,c  =  un  1,  and  once  with  a  =  (pn+1,b  =  (f)n,c  =  (f)n  1. 
This  gives 


^(V+1|2  +  \2un+1  -un\2^j  -  ^(\un\2  +  \2un  -  u"-1^ 

+  -^-\un+1-2un+un~1\2 
4  At  1 

+  i^(V+1l2  +  |20n+1  - ^l2)  ~Ai(m2  +  |20n  ■  1|2) 


+  —  Un+1  -2d)n +  (f)n~1\2 
4A  W  v  1 

+  K+1\2Al  +  \<t>n+1\2A2  +  (C(20n  -  0n-1),wn+1)  -  (CT(2un  -  n"-1), 
=  (/n+1,«n+1)  +  ((7"+1^B+1). 


(7) 


The  second  key  is  to  rearrange  the  coupling  terms.  We  use  the  skew-symmetry 
of  the  coupling  term  and  the  polarization  identity  (Lemma  3)  to  write  it  as 
follows: 


Coupling  =  (C(20n  -  (/U"1),  un+1)  -  ( CT (2un  -  Z"1),  0n+1)  (8) 

=  -(C'(0n+1  -  2 cj)n  +  0n_1),  un+1)  +  (CT(un+1  -  2un  +  u"-1),  (f)n+l) 

=  -W|r+1-20“+r“T  -  A(|„”+ii^cx 

-  j-^|u”+1-2M"+ti”-I|2  -  At|0”+1||rc  +  Kn+I. 

Then  (7)  and  (8)  give 

^(|m”+T  +  |2m"+1  -u”|2l  -  ^(|m”|2  +  |2u"-u”-‘|2') 

+  ^(ir+T  + 12^."+1  -  r  r)  -  ^(in2  +  |20n  -  f-‘i2) 

+  K+‘fi,  +  \r+1\%  -  A«| U"+%CT  ~  A t\r+1\2Crc  +  K"+1 
=  { fn+1,un+I )  +  (gn+1, 4>n+1)- 


Using  again  the  polarization  identity  yields 
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which  by  summation  implies  the  stability  result 
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4  Numerical  verification  of  the  Theorems 


We  give  two  numerical  tests  that  confirm  the  theory  (showing  in  particu¬ 
lar  that  the  restriction  (1)  is  sharp).  The  examples  also  illustrate  that  there 
are  cases  where  each  method’s  time  step  restriction  is  better  than  the  other 
method. 

In  all  test  cases,  the  initial  conditions  are 


and  w1,  (p1  are  computed  using  the  implicit  backward  Euler.  We  take  f  —  g  —  0, 
so  that  any  growth  in  the  energy  is  an  instability. 


Test  1.  In  the  first  case  the  matrices  are 


30  0  \  2  3 

+  =  I  ,  C  = 

0  50  I  l  4  5 


yielding  the  following  time  step  restrictions 

AtcNLF  =  0.1361,  AIbdfab  =  0.2990. 


With  the  time  step 

At  =  0.99  *  AIcnlf 

both  methods  are  observed  to  be  stable,  Figure  1).  With  the  time  step  At  = 
1.01*  A£cnlf  the  CNLF  approximations  exhibit  growth  and  thus  are  unstable 
Since  1.01  *  AtCNLF  <  At bdfab  the  theory  predicts  BDF2-AB2  to  be  stable 
and  this  is  indeed  seen  in  Figure  2. 
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Fig.  1.  Both  methods  stable,  as  predicted. 
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Fig.  2.  CNLF  unstable,  BDF2-AB2  stable,  as  predicted. 


Test  2.  With  matrices 


Ai  — 


Ao  — 


C  = 


the  time  step  restrictions  are 

AtcNLF  =  0.1361,  AtBDFAB  =  0.0299. 


With  time  step  At  =  .99*  AtcNLF  the  CNLF  converges,  while  with  BDF2-AB2 
the  solution  is  unstable,  Figure  3. 
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Fig.  3.  CNLF  stable,  BDF2-AB2  unstable,  as  predicted. 
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